#include "nutation.h"

// IAU2000B 章动序列
const double nuTab[] = {
    0, 0, 0, 0, 1, -172064161, -174666, 33386, 92052331, 9086, 15377,
    0, 0, 2, -2, 2, -13170906, -1675, -13696, 5730336, -3015, -4587,
    0, 0, 2, 0, 2, -2276413, -234, 2796, 978459, -485, 1374,
    0, 0, 0, 0, 2, 2074554, 207, -698, -897492, 470, -291,
    0, 1, 0, 0, 0, 1475877, -3633, 11817, 73871, -184, -1924,
    0, 1, 2, -2, 2, -516821, 1226, -524, 224386, -677, -174,
    1, 0, 0, 0, 0, 711159, 73, -872, -6750, 0, 358,
    0, 0, 2, 0, 1, -387298, -367, 380, 200728, 18, 318,
    1, 0, 2, 0, 2, -301461, -36, 816, 129025, -63, 367,
    0, -1, 2, -2, 2, 215829, -494, 111, -95929, 299, 132,
    // 其他数据省略...
};

// 中精度章动计算表
const double nutB[] = {
    2.1824, -33.75705, 36e-6, -1720, 920,
    3.5069, 1256.66393, 11e-6, -132, 57,
    1.3375, 16799.4182, -51e-6, -23, 10,
    4.3649, -67.5141, 72e-6, 21, -9,
    0.04, -628.302, 0, -14, 0,
    2.36, 8328.691, 0, 7, 0,
    3.46, 1884.966, 0, -5, 2,
    5.44, 16833.175, 0, -4, 2,
    3.69, 25128.110, 0, -3, 0,
    3.55, 628.362, 0, 2, 0
};

// 章动计算
double *nutation(double t, double zq) {
    static double result[2];
    double t2 = t * t, t3 = t2 * t, t4 = t3 * t;
    double l = 485868.249036 + 1717915923.2178 * t + 31.8792 * t2 + 0.051635 * t3 - 0.00024470 * t4;
    double l1 = 1287104.79305 + 129596581.0481 * t - 0.5532 * t2 + 0.000136 * t3 - 0.00001149 * t4;
    double F = 335779.526232 + 1739527262.8478 * t - 12.7512 * t2 - 0.001037 * t3 + 0.00000417 * t4;
    double D = 1072260.70369 + 1602961601.2090 * t - 6.3706 * t2 + 0.006593 * t3 - 0.00003169 * t4;
    double Om = 450160.398036 - 6962890.5431 * t + 7.4722 * t2 + 0.007702 * t3 - 0.00005939 * t4;
    double dL = 0, dE = 0, c, q;
    int i;
    for (i = 0; i < 77 * 11; i += 11) {
        c = (nuTab[i] * l + nuTab[i + 1] * l1 + nuTab[i + 2] * F + nuTab[i + 3] * D + nuTab[i + 4] * Om) / rad;
        if (zq) {
            q = 36526 * 2 * M_PI * rad / (1717915923.2178 * nuTab[i] + 129596581.0481 * nuTab[i + 1] + 1739527262.8478 * nuTab[i + 2] + 1602961601.2090 * nuTab[i + 3] + 6962890.5431 * nuTab[i + 4]);
            if (q < zq) continue;
        }
        dL += (nuTab[i + 5] + nuTab[i + 6] * t) * sin(c) + nuTab[i + 7] * cos(c);
        dE += (nuTab[i + 8] + nuTab[i + 9] * t) * cos(c) + nuTab[i + 10] * sin(c);
    }
    result[0] = dL / 10000000 / rad;
    result[1] = dE / 10000000 / rad;
    return result;
}

// 计算赤经章动及赤纬章动
double *CDnutation(const double z[3], double E, double dL, double dE) {
    static double result[3];
    result[0] = z[0] + (cos(E) + sin(E) * sin(z[0]) * tan(z[1])) * dL - cos(z[0]) * tan(z[1]) * dE;
    result[1] = z[1] + sin(E) * cos(z[0]) * dL + sin(z[0]) * dE;
    result[0] = fmod(result[0], 2 * M_PI);
    result[2] = z[2];
    return result;
}

// 中精度章动计算
double *nutation2(double t) {
    static double result[2];
    double t2 = t * t, dL = 0, dE = 0, c, a;
    int i;
    for (i = 0; i < sizeof(nutB) / sizeof(nutB[0]); i += 5) {
        c = nutB[i] + nutB[i + 1] * t + nutB[i + 2] * t2;
        if (i == 0) a = -1.742 * t; else a = 0;
        dL += (nutB[i + 3] + a) * sin(c);
        dE += nutB[i + 4] * cos(c);
    }
    result[0] = dL / 100 / rad;
    result[1] = dE / 100 / rad;
    return result;
}

// 只计算黄经章动
double nutationLon2(double t) {
    double t2 = t * t, dL = 0, a;
    int i;
    for (i = 0; i < sizeof(nutB) / sizeof(nutB[0]); i += 5) {
        if (i == 0) a = -1.742 * t; else a = 0;
        dL += (nutB[i + 3] + a) * sin(nutB[i] + nutB[i + 1] * t + nutB[i + 2] * t2);
    }
    return dL / 100 / rad;
}